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Abstract 

Motivated by cell adhesion in hydrodynamic flow, here we study bond formation between a 
spherical Brownian particle in linear shear flow carrying receptors for ligands covering the bound- 
ary wall. We derive the appropriate Langevin equation which includes multiplicative noise due to 
position-dependent mobility functions resulting from the Stokes equation. We present a numer- 
ical scheme which allows to simulate it with high accuracy for all model parameters, including 
shear rate and three parameters describing receptor geometry (distance, size and height of the 
receptor patches). In the case of homogeneous coating, the mean first passage time problem can 
be solved exactly. In the case of position-resolved receptor-ligand binding, we identify different 
scaling regimes and discuss their biological relevance. 
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I. INTRODUCTION 



One of the hallmarks of biological systems is their tremendous specificity in binding 
reactions between receptors and ligands. On the molecular level, a prominent example 
is antigen- antibody recognition, which allows our immune system to react to pathogens 
in a highly specific way. Although traditionally much attention has been devoted to the 
biochemical aspects of receptor-ligand binding, physical concepts are equally important in 
this context. In particular, a physical transport process is required to bring receptor and 
ligand to sufficient proximity for binding. A helpful concept is the notion that transport 
has to lead to the formation of an encounter complex, which then can react to form the 
final receptor-ligand complex-i^ii^ii^. In the language of stochastic dynamics, the formation 
of the encounter complex is a first passage problem which can be treated with appropriate 
tools from statistical physics. In many situations, the transport process is simple diffusion. 
However, more complex situations also exist, like the setup in affinity chips, where ligands 
are transported by hydrodynamic flow into a reaction chamber loaded with receptors^. 

In cell adhesion, the physical transport processes required for specific bond formation 
tend to be even more complex, because here receptors and ligands are attached to surfaces 
and their movement is determined by the dynamics of the objects they are attached to. One 
important example in this context are white bloods cells, which circulate the body with the 
blood flow and whose receptor-mediated binding to ligand-coated walls is usually studied 
experimentally in flow chambera^i^i^i^. In order to fight pathogens in the surrounding tissue, 
white blood cells have to extravasate from the blood vessels. Initial binding is provided 
by transmembrane receptors from the selectin family binding to carbohydrate ligands on 
the vessel walls. Here, the probability to form an encounter complex is determined by the 
translational and rotational movement of the cell as determined by hydrodynamic, thermal 
and other external forces. Similar situations also arise in microbiology, when bacteria adhere 
to the intestinal wall^, in malaria infection, when infected red blood cells adhere to the vessel 
wallsr^^'^, in the initial stages of pregnancy, when the developing embryo adheres to the 
uterus^, and in biotechnology, e.g., when sorting cells on microfiuidic chips^^. 

In this paper, we address this situation theoretically by combining methods from hydro- 
dynamics and stochastic dynamics. In Fig. [1] we show the situation which is theoretically 
analyzed in the following. A spherical particle with radius R moves with hydrodynamic flow 
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in positive x-direction at a height z above a wall. The simplest possible flow pattern is linear 
shear flow with shear rate 7. For the usual dimensions in flow chamber experiments with 
white blood cells, this is the relevant flow proflle. In the absence of external forces, there is 
no reason for the particle to drift towards the wall and the formation of an encounter com- 
plex has to rely completely on thermal diffusion. In many situations of interest, however, 
there are forces pushing the particle towards the wall, e.g., gravitational or electric forces. In 
physiological blood flow, cell density is high and the driving force for encounter is provided 
mainly be hydrodynamic or contact interactions with other cells. For the sake of compu- 
tational simplicity and for conceptual clarity, here we consider the simplest case of a force 
driving the particle onto the wall, namely a constant gravitational force directed in negative 
z-direction. Therefore, we introduce a mass density difference Ap between the particle and 
the surrounding fluid. Again this is the relevant situation in flow chamber experiments, 
which are usually done with a diluted solution of cells, thus ruling out a dominant role for 
cell-cell interactions. Receptors are modeled as patches on the particle surface, while ligands 
are modeled as patches on the boundary wall. The formation of an encounter complex is 
then identifled with the flrst approach of any pair of receptor and ligand patches which is 
smaller than a prescribed capture radius tq. The underlying stochastic process is rather 
complex due to position-dependent mobilities resulting from the hydrodynamic equations. 

In order to solve the corresponding mean flrst passage time problem, here we use com- 
puter simulations of the appropriate Langevin equation. A short report of some of our 
results has been given before^. We start in Sec. [Ill by introducing the relevant concepts 
from hydrodynamics at small Reynolds numbers, in particular the friction and mobility ma- 
trices resulting from the Stokes equation for a rigid particle in linear shear flow above a 
wall. In Sec. IIIII we combine these results with concepts from stochastic dynamics in or- 
der to arrive at a Langevin equation describing particle motion subject to hydrodynamic, 
gravitational and thermal forces. Due to the position-dependent mobility functions, we deal 
with multiplicative noise, that is special care is needed to derive and interpret the noise 
terms. In Sec. IIVI our numerical scheme is applied to a sphere falling in shear flow. The 
comparison of the measured stationary height distribution function with the exact solution 
provides a favorable test for our numerical treatment. In Sec. |V]we show that for the case 
of homogeneous coverage of sphere and wall the mean flrst passage time to contact can be 
solved exactly, again in excellent agreement with our numerical procedure. In Sec. |Vl] we 
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explain why the choice for the initial height is not essential. In the next two sections, we 
present and explain our simulations results, first for movement restricted to two dimensions 
in Sec. IVIII and then for the full three dimensional case in Sec. IVIIII We finally conclude in 
Sec. IIXI by discussing the biological and biotechnological relevance of our results. 



II. FRICTION AND MOBILITY MATRICES 



Due to their small sizes, the hydrodynamics of cells is in the low Reynolds number regime. 
Using a typical cell size L = 10 fim and a typical velocity v = mm/s (that is the flow velocity 
at a distance L to a wall with linear shear flow of rate 7 = 100 Hz), the Reynolds number is 
Re = pvL/f] = 10^^, where p = g/cm^ and rj = 10^^ Pa s are density and viscosity of water, 
respectively. Therefore, we essentially deal with the Stokes equation for incompressible 
fluids: 

r]Au{r) -VP{r) = -T{r), V • u(r) = 0, (1) 

where u(r) is the fluid velocity field, -P(r) is the pressure field and J-'{r) is the force density 
on the fluid by the particle. Here, we use the induced force picture, i.e., the fluid equations 
of motion Eq. ([T]) are extended to the interior of the particle and the particle is replaced by 
an appropriate force density J^(r) acting on the fluidr^. The unperturbed flow field has to 
satisfy the homogeneous version of Eq. ([1]) as well as no-slip boundary conditions at the wall. 
In this paper, we use the simplest possible example, namely linear shear flow, u°° = 'jze^. 
The effective flow field in the region occupied by the rigid sphere reads 

u(r) = (U + ^^ X (r - R)) 0(i? - ||r - R||) , (2) 

where U, Q are the translational and rotational velocities of the sphere, respectively. R is 
the position of its center, R the sphere radius and G the theta step-function. The particle 
is subject to forces and torques which follow from the force density as 

= j J^(r)dr, = J{r-R)x J^(r)dr . (3) 

Because we consider a rigid object, higher moments of the force density are not required in 
our context. For the unperturbed flow at the mid-point of the sphere, we make the following 
definitions: 

1 1 

, (4) 



U°° = u°°(R), ri°° = ^ V X u° 
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r=R 
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where the vector $7°° is called vorticity and the tensor E°° rate of strain tensor. Because we 
restrict ourselves to linear shear flow, all higher moments of the unperturbed flow vanish. 

Due to the linearity of the Stokes equation, a linear relationship exists between the force 
density J^(r) and the driving flow, which is the difference between real and unperturbed 
flows^^. Specifled for the flrst moments of the force density, it leads to the relation 



(5) 



where the shear force = Re : E°° with A : B = tr AB""^. It results from the perturbation 
of the flow by the presence of the wall and vanishes for free flow. The two matrices Ru and 
Re are conveniently written as 
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R-= \ , Re := _ , (6) 



where the ( are the symmetric friction matrices and the superscripts t, r and d stand for 
translational, rotational and dipolar, respectively. In order to obtain the translational and 
rotational velocities of the sphere as a function of the hydrodynamic forces and torques, we 
have to invert Eq. ([5]): 






(7) 



The symmetric matrix M = Ru ^ is called mobility matrix. It is convenient to deflne the 
mobility tensors through 

M ^ R„- J V R-R.^""). (8) 

\i^'' i^" I V/^^'v 

In order to calculate the friction and mobility tensors for the special case of a sphere in 
linear shear flow above a wall, we follow the procedure from Ref.— . The friction tensors C, 
introduced in Eq. (jS]) and the mobility tensors /i introduced in Eq. ([HD are expressed in terms 
of scalar functions together with irreducible tensors formed form the Kronecker symbol 
the Levi-Civita symbol eijk and the normal vector k = e^. The scalar friction and mobility 
functions are not known analytically, but can be obtained to high accuracy by the following 
numerical scheme. One introduces the variable t = R/z, where R is the radius of the sphere 



and z is its height above the wall. Thus, t can take values from the interval [0, 1]. In the limit 
t ^ 0, that is far away from the wall, one can expand the friction functions in powers of t. In 
the limit t — > 1, that is close to the wall, analytical results can be obtained with lubrication 
theory. In order to cover the whole interval, the two limit solutions are matched using a 
Pade summation scheme^^i. More details of this implementation are given in appendix El 
In contrast to the tabulated finite element results from Ref.— , this implementation gives 
correct results for any possible configuration. 

III. LANGEVIN EQUATION 

The motion of a particle subject to thermal, hydrodynamic and direct external forces 
like gravity is called Stokesian Dynamics"^^. In this section we derive the corresponding 
stochastic differential equation {Langevin equation). The Langevin equation will allow us to 
base our statistical treatment on the repeated simulation of individual trajectories. Because 
we are interested in the over-damped (Stokes) limit, we can neglect inertia in Newton's 
second law: 



where — F , F and F are hydrodynamic, direct and thermal forces acting on the sphere. 
An analogous balance exists for the torques. For the following, forces and torques as de- 
scribed above are united in one symbol. For example, from now on the symbol F denotes 
(F,T), a six-dimensional vector comprising force F and torque T, and U denotes the six- 
dimensional particle translational/rotational velocity vector. 

In the absence of Brownian forces, F^ = and F'^ = F"^. Inserting this into Eq. ([7]) then 
gives 



and the particle trajectory can be found with a simple Euler algorithm as X(t + At) = 



In the presence of Brownian motion, the situation is more complex, because thermal noise 
leads to terms of the order At^/^ and special care has to be taken to include all terms up to 
order At. Due to the fluctuation-dissipation theorem, for our problem Gaussian white noise 



_-pH 4. + F^ = 



(9) 




(10) 



X(t) + UAt + C(At2). 
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reads 

iSt) = 0, (gtgf) = 2kBTM6{t - 1') . (11) 

Here, the subscript t corresponds to the fact that the thermal force g is a random process. 
The left part of Eq. (ITT!) states that the forces that the fluid exerts on the particles are 
equally distributed in all directions so there is no net drift due to thermal fluctuations. The 
right part of Eq. (fTTj) states that forces at different times are not correlated, which is a good 
approximation because the diffusive forces act on a much faster time scale than the hydrody- 
namic forces. Because the mobility matrix M is position-dependent, we deal with so-called 
multiplicative noise. Since the 5-correlation in Eq. flTTl) can be considered to be the limit 
of a process with an intrinsic time scale for thermal relaxation, which is much faster than 
the time scale of hydrodynamic movement, the Stratonovich interpretation of the stochas- 
tic process is appropriate^^. This means that for each time step, the mobility functions 
have to be evaluated at X(t + (1/2) At) (rather than at X(t) as in the Ito interpretation). 
The Stratonovich interpretation also implies that the rules for integration and coordinate 
transformation are the same as for the Riemann integral in non-stochastic calculus. 

The presence of the thermal noise Eq. (ITT!) converts the position function X(t) into a 
random process X^. Multiplicative noise can result in additional drift terms. We therefore 
write the Langevin equation as 

9tXt = U°° + M(F^ + F^) + A;BrY + gf, (12) 

where in comparison to the deterministic equation Eq. ffTOl) we have added both the Gaussian 
white noise gf (to be interpreted in the Stratonovich sense) and some drift term Y. The 
drift term Y can be derived by requiring Eq. f|T2|) to be equivalent to the appropriate 
Smoluchowski equation. The details of these calculations are given in appendix [Bl The 
result is 

Y=BVB'^, M = BB^, Yi = BikidiBik), M,, = B.^B.^ . (13) 

For additive noise, that is for position-independent mobility functions, the additional drift 
term would vanish. In the case of position-dependent mobility matrices, the noise term 
gf alone would lead to a drift of the particle towards regions of lower mobility (that is 
towards the wall, where mobility vanishes due to the no-slip boundary condition). This 
drift, however, is exactly compensated by the additional term Y. 
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For the following, it is useful to non-dimensionalize Eq. (fT2|) . For length, the natural scale 
is sphere radius R. For time, we use QTtrjB? /ksT , which is the time needed to diffuse the 
distance R. For force, we use Gii'qR'^'y, the Stokes force at velocity -R7, that is in linear shear 
flow a distance R away from the wall. The scalar friction and mobility functions appearing 
in M, Re and Ru, also become dimensionless as explained in appendix [XI The Langevin 
equation Eq. (fT2|) now reads 

9tXt = Pe (U°° + M(/F^ + F^)) + BVB^ + gf , (14) 

where the Peclet number Pe = GnriR'^'y /ksT measures the relative importance of determin- 
istic to Brownian motion. In the limit Pe the particle only exhibits diffusive motion 
and in the limit Pe ^ 00 it is no longer subjected to diffusion. The second dimensionless pa- 
rameter / = ||F'^||/67r?7i?^7 measures the relative importance of direct forces/torques versus 
the shear force/torque. Measuring the time in units of the diffusive time scale is appropriate 
for Peclet numbers of order ten or less. For simulations with larger Peclet numbers it is 
more suitable to scale time with the inverse shear rate 7"^. This has the effect of dividing 
Eq. ([H]) by Pe. 

In order to solve Eq. (fT^ numerically, it has to be discretized with respect to time. The 
appropriate Euler algorithm can be derived by first rewriting Eq. (HM in the Ito-version, 
which adds another drift term to the equation. As explained in appendix [Ul the two drift 
terms together lead to the result 

dt^S. = Pe (U°° + M(/F^ + F^)) + VM + g,^ . (15) 

Its discretized version is simply 

AX = [Pe (U°° + M(/F^ + F^)) |^ + VM|,] At + g(At) + 0{At^) . (16) 

This final result has been derived before in a different way by Brady and Bossis^i. For 
vanishing shear flow, it also agrees with the classical result by Ermak and McCammon^^. 
In appendix O) we describe the algorithms used to implement Eq. (fT6|) , in particular the 
algorithm to generate the thermal forces g(At). 

IV. SPHERE FALLING IN SHEAR FLOW 

As explained in the introduction, we consider a sphere whose density is slightly larger 
than that of the fluid. Due to this density difference Ap a constant drift towards the wall 
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exists. As we will see later, this drift ensures that on average the sphere will bind to the 
wall in finite time. The two independent parameters defined in Eq. f|T^ for this model 
system are Pe and / = {2RApg)/{9T]'y), with the earth acceleration constant g = 9.81 
m/s^. For later considerations, it is convenient to introduce also the parameter Pe^ = f Pe, 
which we call the Peclet number in z -direction. Pe and Pe^ represent the strengths of the 
hydrodynamic and gravitational forces in respect to the thermal force, respectively. Out of 
the three parameters Pe, f and Pe^, only two are independent, because / = Pez/Pe. 

We first consider the path of a sphere falling in shear fiow after it has been dropped 
at some initial height. Fig. [2] illustrates the effect of the Peclet number by showing some 
representative simulation trajectories. For Pe = oo the motion of the sphere is purely 
deterministic and only governed by the parameter /. In the diffusive limit Pe = 0, the 
sphere makes a pure random walk (except for the drift in 2;-direction due to the gravitational 
force) . 

As the mobility matrix does only depend on the height of the sphere above the wall (cf. 
appendix El), the motion in the z-direction is independent of the position in the (x,?/)-plane 
and the orientation of the sphere. Therefore, it can be treated separately. The probability 
density "^{z, t) for the sphere to be at height z at time t is the solution to a one dimensional 
Smoluchowski equation 



This equation cannot be solved analytically as the mobility function Mzz is not known in 
closed form. In Fig. [3] we show numerical solutions obtained by simulating the equivalent 
Langevin equation. One clearly sees that first the 5-function at t = is broadened due 
to diffusion and then develops into a stationary solution which has its maximum at the 
wall. This stationary solution has a simple analytical form which follows from Eq. ( ITTj) by 
integrating = 0: 



Thus, the stationary solution is simply the barometric formula, as it should be for ther- 
modynamic reasons. We also find that the first two moments (mean and variance) are the 
same: 



dt^{z,t) = -dzJz, Jz = -Mzz{dz^ + Pez^). 



(17) 



^siz) = Pe,e-^^^(^-i) . 



(18) 



{z-i) = y/{z^)-{zy 



1 



(19) 



Pe. 
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In the limit of vanishing gravitational force (Pe^ 0), the probability distribution becomes 
flat and the probability of finding the sphere does not peak at the wall anymore. 

V. FIRST CONTACT WITH HOMOGENEOUS COVERAGE 

If the sphere and the wall are homogeneously covered with receptors and ligands, respec- 
tively, an encounter complex is established whenever the sphere comes sufficiently close to 
the wall. The mean time which elapses after the sphere is set free at some initial position 
until an encounter complex is established is then identical with the mean first passage time 
(MFPT) for a sphere dropped at initial height zo to reach the height zi. Note again that 
the motion in z-direction is independent of the values of the other coordinates. For a par- 
ticle diffusing in an interval [zi, b], with zi being an absorbing boundary and b a reflective 
boundary, the MFPT T to reach Zi when started at 2; G [zi, b] is the solution to the following 
ordinary differential equation^^ 

A{z)d,T{z\zi) + D{z)dlT{z\zi) = -1, T{zi\zi) = 0, d,T{z\zi\^, = 0. (20) 

In our case, b = 00. The drift term is A{z) = —Pez<y'^^{l/z) + (9^q;**(1/z) and the diffusive 
term D{z) = M^z = «**(l/-2), where q;**(1/z) is a scalar mobility function as explained in 
appendix [Al The general solution to Eq. (!20l) 

zo / 00 \ / z \ 

^(^°|^') - [ '^W) [I ''mj • ^ [I'^m) ■ '''' 

This, can be reduced up to an integral over d**(l/z): 

zo 

r(^„k,) = ^/d.^. (22) 

Zl 

Thus the dependence of T{zo\zi) on Pe^, the only parameter in this problem, is obtained 
exactly. It is important to note that the compact form for the MFPT in Eq. (122|) is a result 
of the constant vertical force. For a more general vertical potential force F± = —dzV{z) 
with a potential V, Eq. (1211) can be reduced to 



T{zo\zi) = 
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This equation shows that the potential must satisfy the condition limy^ooiV (z) — V{y)) — >■ 
— oo for the MFPT to be finite. This holds true, e.g., for the gravitational force studied here 
or for the interaction of a charged object with an oppositely charged wall, but not, e.g., for 
a Lennard- Jones potential. 

The integral Eq. (1221) over the scalar mobility function can easily be calculated nu- 
merically as d** behaves well in the full range of z. In fact a^^{t) can be approximated by 
its leading term from the lubrication analysis, i.e., ^ 1 — t. We then find 



A numerical analysis shows that the approximation Eq. fl24l) deviates only by a few percent 
from the exact solution Eq. (!22|) . Thus, T(2;o|-2;i) is logarithmically divergent if the absorbing 
point is close to the wall, zi ^ 1, and linearly divergent if the starting point is at infinite 
height, zo ^ oo. 

For a sphere homogeneously covered with receptors each having a capture radius ro, 
the mean time for forming an encounter complex is T{zo\l + tq). This time will serve as 
a useful limiting result in some of the considerations presented in the next sections. The 
exactly known result Eq. fl22|) provides also a good test for the algorithm we implemented. In 
Fig. the MFPT obtained from simulation experiments and from quadrature of Eq. ( |22l) are 
compared. The two results agree very well (see appendix [D] for a discussion of the statistical 
and systematic errors of the simulation results). In Fig. IDd we show the numerically obtained 
distribution of first passage times. One clearly sees that the larger Pe^, the stronger they 
peak around the mean. 

We conclude the case of homogenous coverage by noting that in order to obtain dimen- 
sionalized results, one has to multiply the MFPT by the diffusive time scale GnriR^/kBT. 
This result does not depend on shear rate 7 because vertical and horizontal motion are 
decoupled and rotational motion is not relevant here. However, it depends on viscosity 77, 
which sets the time scale for vertical motion. If one switched off thermal fluctuations, the 
falling time would be exactly the same as the MFPT from Eq. (!22|) . but this is a special 
result for constant force and not true in general. If one removed the wall, the translational 
symmetry in z-direction would not be broken and the MFPT would be T = (zq — zi)/Pez, 
that is the logarithmic term in Eq. ( 12^ would be missing. 




(24) 
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VI. EFFECT OF INITIAL HEIGHT 



We now turn to spatially resolved receptor coverage, that is we consider a sphere which 
is covered by Nr equidistantly spaced receptor patches. For the moment being, the wall is 
still considered to be homogeneously covered with ligands. The MFPT T{6,x\C) now will 
depend on the initial position x = {x, y, zq) and the initial orientation 6 as well as on the 
absorbing boundary C in diffusion space. The latter is given by the special receptor and 
ligand geometry. In an experimental setup with linear shear flow it is possible to measure 
only particles which have been initially at a certain height. This is due to the fact that 
their average velocity as obtained from the solution of the Stokes equation Eq. ([7]) depends 
on their height in a unique way*. However, it is almost impossible to prepare a certain 
initial orientation 6 or {x, ?/)-position relative to the ligands. Therefore, the quantity of 
interest to us will be a MFPT which is averaged over all possible initial orientations 6 
and all initial positions (x, y), which will be denoted as {T{9, x\C))ff^^ The dependence of 
(T(x, 0\C))f^^^ on the initial height for zq > I+vq can be derived exactly. For homogeneous 
ligand coverage the quantity of interest is 

Je 

where C is the absorbing hyper-surface in (6', 2;)-space and Vg a normalization constant. 
Absorption is only possible if 2; < 1 + ro, thus if we look at some intermediate height 

-2^0 > -2^m > 1 + '"o; then 

T{lz^\C) = T0,z,\z^) + j d^emP{Om\0)T{L,Zm\C), (25) 

where p{6m\(^m) is the conditional probability to pass the height Zm with the orientation 6m 
when starting with the initial orientation 6 at zq. T{6, zo\zm) is independent of the initial 
orientation and can be calculated by means of Eq. fl22l) . Now averaging Eq. fl25l) over the 
initial orientation gives 

T{L,Zm\C) (26) 

= T{Zo\Zm) + r^ I (f6mT{em,Zm\C)=T{zo\Zm) + {T{6m,Zm\C))g^. 

Thus, if the orientation-averaged MFPT is known for some initial height Zq > 1 + tq, then 
the MFPT for any other initial height z'q > 1 + Vq can be calculated by means of equations 
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{n9,zo\c))g=nzo\zm) + 



Eq. (l26|l and Eq. (122|) . In Fig. [5l this result is verified by simulations for the two-dimensional 
case, that is the sphere can only move in the x — z-plane and rotate only around the y-axis, 
compare Fig. [H Due to the decomposition Eq. (12^ . the initial height is not essential. In 
the following, we therefore will always use the value zq = 2, that is the sphere has to fall by 
one radius until it hits the substrate for the first time. 

VII. MOVEMENT IN TWO DIMENSIONS 

We now study the effect of shear rate for heterogeneous receptor distribution if the sphere 
is restricted to move only in two dimensions. Then, the receptor patches can be equidistantly 
distributed over the circumference as illustrated in Fig. [61 Each receptor patch has a capture 
height of To and a width of 2rp. The 2D receptor density is then pr = NrTp/ii. Orientation is 
now represented by a single angle 9. The absorbing boundary C is illustrated in Fig. O For 
each receptor patch, binding can occur over a range 29o, which consists of two parts. The 
inner part is valid already for rp = and reflects the overlap due to a finite tq. The outer part 
is results from a finite r^. Together this leads to 6q{z) = arccos(z/(l + ro)) +rp. The receptor 
patches establish a periodicity with period 6s = 2'K/Nr. As the number of receptor patches 
grows, this period decreases and one finally achieves overlap. Then, encounter becomes 
possible for all values of 6, that is we are back to the case of homogeneous receptor coverage. 
In our case of non-homogeneous coverage, the MFPT depends on Pe, Pcz, Nr, ro, Vp and zq. 
For the following simulations Vp = tq = 10~^, Pe^ = 50 and 2:0 = 2 is chosen unless other 
values are explicitly mentioned. 

Fig. [7^ shows the MFPT as a function of the Peclet number Pe. Note that in the log-log 
plot, an apparent plateau appears at small value of Pe, although in a linear plot there would 
be monotonous decay. Three regimes can be distinguished. For Pe ~ {diffusive limit) the 
transport by the imposed shear flow is negligible and only diffusive transport is present. For 
very large values of Pe, (T)^ plateaus at the value given by Eq. fl22|) independent of A^^,- In 
this limit the time for rotation to any certain orientation is negligible compared to the mean 
time to fall down close to the wall, therefore, the result for rotational symmetry is recovered. 
Between these two limits the MFPT decreases monotonically with increasing Pe. Fig. Wp 
shows the data from Fig. [7^ plotted as a function of the receptor density p^. (x N^. The larger 
Pe the less pronounced is the dependence of (T)^ on N^. For Pe ~ 0, however, {T)q strongly 
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depends on A^^. The latter relation is better illustrated in Fig. [Tt. There, at Pe ~ 0, {T)o 
is shown for a wide range of N^. The simulations were done for fixed patch size Tp but for 
four different values of the capture radius tq (cf. Fig. E]). For pr I, {T)g reaches the value 
given by Eq. fl221) . As described by Eq. fl221) . {T)g is the smaller the larger tq is. An increase 
in the number of receptor patches N,. leads to a strong decrease for the MFPT, however, no 
special scaling behavior can be observed. It is remarkable that the limiting value for the case 
of homogeneous receptor coverage is already reached for pr ~ 10~^. The larger the capture 
radius tq the more pronounced is this effect. This can be understood by observing that the 
effective patch size as given by the angle 60 > Vp (see Fig. [6]) is monotonically increasing 
with increasing tq. 

We next try to qualitatively understand the effect of shear rate for the simulation results 
shown in Fig. [7^. In general, it is very hard to separate the effects of diffusion and convection. 
The time for binding at Pe ^ is determined purely by diffusion effects and will be denoted 
by Td. As shear flow increases, the rotation of the sphere is increasingly dominated by 
convection. We now derive a convection time Tp which competes with the diffusion time 
Td at large Peclet number. For very large Peclet number, we expect the MFPT to be the 
sum of the homogeneous result from Eq. fl24|) plus this additional time Tp. An important 
question then is at which Pe the convection time Tp become smaller than the diffusion time 
Td. 

On order to estimate Tp, we note that the main effect of increased shear rate is faster 
rotation in the direction of flow. Once a receptor has rotated by an angle Og = 271 /Nr such 
that it opposes a ligand on the substrate, there is some probability p that the sphere is at the 
correct height that an encounter can occur. If no encounter occurs with the complementary 
probability 1 — p, the sphere has to rotate about another angle 6s until the next receptor 
points downwards. Supposing the time, 2to, to rotate about the angle 9s is large enough that 
there is no correlation between the height of the sphere before and after the rotation, then, 
an encounter occurs again with probability p (therefore this analysis also does not hold at 
very large Pe). Thus, the mean time Tp for encounter is 

Tp = pto + (1 -p)(p3to + {l-p){p5to + (1 -p)(. . .))) 

= pto 5^(2z + !)(!- pY = to'-^ ^ ^ , (27) 

U p p 

where the series has been summed up by means of the geometric formula. In the last term 
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we assumed that the probabihty p for the proper height is small due to a small capture 
distance tq. It follows from the stationary probability distribution "^siz) given by Eq. f[T8l) : 

"l+ro 



P 



/i+ro 
dz^siz) = 1 - 6"^'^^'^° ^ Pe^ro . (28) 



The time to to rotate about half of the angle 6s is approximately to = 6s/ Pe. Therefore, we 
get 

„ 47r 



(29) 



NrPePe^ro 

In this analysis, the convection time Tp scales inversely with the number of receptor patches 
Nr and the Peclet number Pe. As Pe increases, Tp gets smaller than To and then dominates 
the overall outcome. Comparing Eq. (l29l) to the simulation data for Pe ~ shows that this 
crossover occurs in the range Pe ~ 10^ — 10^ and that the corresponding value of Pe increases 
with increasing receptor number Nr, exactly as observed in the simulation data over the full 
range of Pe. However, the exact scaling of this data is not ~ l/A^r for large Pe as predicted 
by Eq. fl29l) . In practice, the decay is somehow slower due to correlations between the height 
of the sphere at two successive instances of a receptor pointing downwards, which we have 
neglected in our analysis. 

We briefly comment on the effect of the downward driving force, that is Pe^. Above, we 
have found that in two cases, homogeneous coverage from Eq. (122!) and convection-dominated 
rotation from Eq. fl29|) . the MFPT scales inversely with Pez- This scaling behavior is indeed 
found in the simulations, except that for very large values of Pez, the MFPT approximates 
a constant value (data not shown). The reason is that the larger Pez, the smaller the mean 
time to fall below the height z = 1 + vq. As indicated by Eq. ( JTSll . then the sphere stays 
below this height until an encounter occurs. This implies that in this limit, the MFPT 
depends only on rotational motion and the falling motion is irrelevant. 

We now introduce spatially resolved ligands into the 2D-model. Fig. [8^ shows the model 
definition: the ligand patches are considered to have the same radius = Vp as the re- 
ceptor patches and they are located at a distance d from each other. This results in a 
one-dimensional ligand density given by pi = 2rd/d. The mean first passage time will now 
also depend on the initial x-position, T = T{zo,9,x\C), where C is the hypersurface in 
{z, 9, x) space where a receptor patch touches a ligand patch. But similarly as in the above 
section in regard to initial orientation, the dependence on the initial x-position is of mi- 
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nor interest and therefore, we will discuss the MFPT averaged over the initial position and 
orientation, denoted by {T)e,x- 

Fig. shows that by varying the Peclet number we can identify the same three regimes 
for all ligand-densities as before. For Pe ^ in the limit of pure diffusive transport, {T)q^x 
approaches a finite value, depending on and pi. With increasing Pe, (T)^ ^; decreases 
monotonically and finally for Pe — > oo reaches the value of the MFPT in the limit of 
homogeneous receptor and ligand coverage. In contrast to above, however, in this limit 
the shear flow not only restores rotational invariance of the sphere, but in addition also 
translational invariance of the substrate. 

Fig. provides more details for (T)^ .j. as a function of pi in the diffusive limit (Pe ~ 
0). We find that in the range 0.1 < < 1 the MFPT is almost not affected by ligand 
concentration: as long as the ligand patches are sufficiently close to each other, a receptor 
patch touching the wall will most probably find a ligand before diffusing away again. The 
situation changes completely with small ligand density. For pi the averaged mean first 
passage time {T)o^x scales with the ligand density pi as {T)o^x oc 1/pf oc d?. This can be 
understood by calculating the position- averaged MFPT (T)^. for a particle diffusing in an 
interval [0, d] with diffusion constant which gives {T)^ = d^/l2D. This suggests that the 
quadratic scaling with d results from the diffusive motion between adjacent ligand patches. 
Fig-[9t> summarizes our results for the dependence of the 2D MFPT (T)^ .j. on ligand density 
pi and receptor density pr in the diffusive limit. Clearly there exists a large plateau around 
the value for the case of homogeneous coverage pr = pi = 1. This implies that if ligands and 
receptors patches are not too strongly diluted, the mean encounter time is still close to the 
optimal value given by Eq. (122!) . On the other hand if the number of receptor and / or ligand 
patches is highly reduced the mean encounter time is strongly increased. 

VIII. MOVEMENT IN THREE DIMENSIONS 

We finally turn to the full 3D-situation, that is the sphere may diffuse about all three 
described by Eq. f|T6|) and Eq. flC4|) . Receptors are located in spherical patches 
which are randomly distributed over the sphere. Each receptor patch has a radius Vp and 
a height (capture length) tq. That is the appropriate generalization of the situation shown 
in Fig. [6] for the 2D-case. Thus, for A^^ receptor patches the receptor density is pr = 
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27rA^r(l — cos(rp))/47r ^ NrVp/A (for ^ 1). In contrast to the preceeding sections where 
the receptor patches could be regularly distributed over the circumference, this is no longer 
possible on the surface of a sphere. Therefore, we distribute the patches randomly over the 
sphere with equal probability for each position, with a hard disk overlap algorithm making 
sure that no two patches overlap^. One has to bear in mind that then for small A^^ two 
different distributions may have slightly different binding properties. This effect becomes 
weaker for larger Nj., therefore in the following we will only use A^^ > 10. The quantity 
we measure in our simulations is now (T)^ in the case of homogeneous ligand coverage and 
('^)e{xy) ^^^^ ^^^^ ligands are located in spherical patches on a 2D-lattice. Thus, 

we average the MFPT over the initial orientations and positions as explained above. 

In order to explore the dependence of (T) ^ on A^^ and Pe we first simulated the receptor 
ligand encounter in the case of homogeneous ligand coverage pi = 1. In order to average over 
the initial positions we started each run with a randomly chosen initial orientation. After 
100 runs we generated a new distribution, thus averaging out also the effect of different 
receptor distributions. In order to achieve reasonable statistics, we typically used 100,000 
runs. Our results are shown in Fig. [TOk . Again we find three different regimes as a function 
of the Peclet number Pe. This proves that qualitatively the basic results of the 2D-treatment 
remain valid in 3D. However, in detail there are important differences. In contrast to the 
2D results presented above, (T)^- in the limit Pe ^ cxo is no longer given by Eq. fl22|) if A^^ 
is small. That is due to the fact that for Pe —>■ oo the receptor patches effectively behave 
as ring-like structures. The rotation of such a ring about the x- or y-axis is not affected by 
Pe and thus still depends on diffusion. For large A^,. the rings cover the whole sphere and 
for Pe ^ oo (T)g' is again given by Eq. (!22|) . 

In Fig. [TOb we plot the Pe — > limit of (T)^' as a function of the number of receptor 
patches A^^, for different values of the capture radius tq. The fitted straight line for vq = 10~^ 
shows that (T)^- approximately behaves like (T)^' oc 1/Nr. Neglecting effects of curvature, 
the average distance between two receptors patches is oc {ATc/NrY^"^ and the mean time to 
diffuse that distance is td oc (x 1/Nr. This provides a simple explanation for the observed 
scaling behavior. For high Nr, the MFPT reaches a plateau value, given by Eq. fl22|) . This 
plateau value depends on tq and is the smaller the larger rg. Also the crossover from the 
asymptotic behavior at small A^^ to the plateau at large A^^ is shifted with increasing capture 
radius tq towards smaller A^^- 
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In Fig. [TOb we show the effect of a finite hgand density pi at Pe ~ 0. For the simulations 
we distributed the hgands in circular patches of radius = 0.01 on a quadratic lattice with 
lattice constant rf, thus, resulting in a ligand density pi = vrr^/rf^. In our implementation, 
the intersection between the receptor patch and the wall is approximated by an appropriate 
circle, because it is easy to check if this circle overlaps with the ligand patch. The fits given 
in Fig. [TOb show that for small pi, the MFPT scales as (7")^ oc 1/pi oc cP. Because the 
curves for different Nr appear to be rather similar, in the inset we plot the ratio of different 
pairs of these curves. As this results in approximately constant plateaus, we conclude that 
the scaling with ligand density is hardly effected by Nr. As in 2D, the inverse scaling with 
ligand density can be understood in simple terms by noting that the MFPT to diffusional 
capture scales like cP. At a coverage around 0.01, saturation occurs as it did for receptor 
coverage. 

We finally discuss the influence of the receptor geometry described by the parameters tq 
and Tp. Because Pe changes the MFPT in a monotonous way, it is sufficient to study the 
diffusive limit Pe ~ 0. Fig. [TTk and b show (T)^- as a function of for ro = 0.001 and 
ro = 0.01, respectively. In order to obtain smooth curves, in this case only one receptor 
distribution was used for all runs. We find that the curves can be fitted well to the function 

(^(^p))e = 7-^ + Tizo = 2\zo = 1 + ro), (30) 

where the second term is the homogeneous result from Eq. ( l22l) . This means that even for 
vanishing receptor size rp — > the MFPT remains finite. This makes sense because above 
we have shown that the effective patch size is determined both by Vp and tq. In detail. Fig. [6] 
showed that capture occurs over the sohd angle 26^0 with 9o{z) = arccos(2;/(l + ro)) + r^. 
For small ro and r^, this allows us to define an effective patch size 

r^-^-^ = arccos((2;)/(l + ro)) + Vp ^ arccos(l - ^ro) + r^ ^ ^/tq + rp, (31) 

where we have used (z) = 1 + ro/2. Suppose now that the sphere diffuses over the time 
until a receptor patch points downwards, then it may encounter a ligand with a probability 
p that is given by the normalized area of one effective receptor patch: 

P= ^(1 -cos(r;^0) ~ liV^ + rpY ^ \v^{lv^ + rp). (32) 

If no encounter occurs, the sphere has to diffuse again a time until the next encounter can 
occur. This leads to the mean encounter time T = td/p- Putting everything together gives 
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Eq. (l30l) with a = 2trf/(yro) and b = ^^/rQ■ If checked against our simulation results, we 
indeed find that the fit parameter b is an increasing function of tq, but varies only slightly 
with Nr- The fit parameter a scales approximately as ~ l/N^ and varies with tq, also 
consistent with the above analysis. In Fig. [TTb (T)^- is plotted as a function of Vp for several 
values of tq and A^^ = 30. One clearly sees that increasing has a much smaller impact on 
(T)g>than a comparable increase in ro, which is qualitatively well described by the preceeding 
analysis. 

In Fig. [TTb and b the receptor density is varied over almost four orders of magnitude by 
changing r^, but the largest measured decrease for (T)g is only by a factor four. In contrast, 
an increase of the receptor density by one order of magnitude due to ten-fold more receptor 
patches leads to a decrease of (T)^by almost also one order of magnitude. However, this is 
only true as long as A^^ is not too large, as for large Nr {T)g saturates at the limiting value 
of homogeneous receptor coverage (cf. Fig. [TOb). The crossover from the behavior to 
the saturation should take place when the average distance between two receptor patches 
d! ~ (47r/A^r)^^^ becomes comparable to the size of one receptor patch. This corresponds 
to Tp^^ ~ {AT^/Ny)^^"^ or Nr ~ 47r/(y¥o + Tp). This estimate predicts that the crossover 
takes place between several tens to several hundreds of receptor-patches, depending on tq, 
in agreement with the data shown in Fig. [TOb . 

IX. SUMMARY AND DISCUSSION 

In this paper we have calculated the mean first passage times (MFPT) for initial en- 
counter between spatially resolved receptors on a Brownian particle in linear shear fiow and 
spatially resolved ligands on the boundary wall. Our main results were obtained by repeated 
simulations of the discretized Langevin equation Eq. fll6p . Each data point shown corre- 
sponds to at least 100,000 simulation runs. It is important to note that these simulations 
are very time consuming because we resolve objects of the size of 10~^i?, that is for /xm-sized 
particles we resolve the nm-scale. 

In general, we found that the MFPT was always monotonically decreased when the Peclet 
number was increased. That means that a particle which is covered with receptors in a way 
that it binds well to ligands already in the diffusive limit is even better suited to initiate 
binding at finite shear rate. In our simulations we modeled the receptor geometry using 
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three parameters: the number of receptor-patches Nr, the radius of the receptor patches 
Tp, and the capture radius tq. The efficiency of binding is mainly increased by Nr, but 
only up to a saturation value of the order of hundred. An increase of leads only to a 
weak enhancement of binding efficiency. The influence of tq to the MFPT is threefold: i) 
it reduces the mean falling time, ii) it increases the effective patch size, and iii) according 
to the stationary probability distribution for the 2;-direction, it becomes more probable for 
the sphere to be within the encounter zone when tq is increasing. An additional but more 
indirect effect of receptor protrusions is that the further the cell is away from the wall, the 
faster it can rotate (even in the diffusive limit) due to the larger mobility. As shown by 
Eq. (1261) rotations play a role only within binding range, i.e., for z < 1 + tq. Therefore, a 
large tq lets the cell also benefit from faster rotations. Summarizing our findings in regard 
to receptor geometry we conclude that the most efficient design for particle capture under 
fiow is to cover the particle with hundreds of receptor patches (A^^ above threshold), each 
with a rather small area (small r^), but formed as a protrusion (large ro). 

Indeed, this strategy seems to be used by white blood cells, which have evolved intriguing 
mechanisms both on the molecular and cellular scale in order to adhere effectively to the 
endothelium under the conditions of hydrodynamic fiow. The typical size of white blood cells 
is i? ~ 5 fim and they are covered with a few hundreds of protrusions {microvilli) with the 
receptors (most notably L-selectin) localized to the microvilli tips2^. In general, the microvilli 
of white blood cells are much more complex than the parameter ro in our model: they are 
rather long (typical length 350 nm, that is R/15) and have their own physical properties 
(e.g., very fiexible in the transverse direction and viscous in the longitudinal direction)^. 
Nevertheless, it is striking that elevation of the receptors above the main cell surface seems 
to be a major design principle for white blood cells. In fact, the same strategy appears to be 
used also by malaria-infected red blood cells, which are known to develop a dense coverage 
with elevated receptor patches (knobs) on the cell surface^»^>^. A typical value for the cell 
radius is 3.5 fiirir^. The knobs have a typical height of 20 nm, a radius of about 90 nm and 
a distance of 200 nm (for red blood cells infected by single parasites)^. This dense and 
elevated coverage suggests that like the white blood cells, the malaria-infected red blood 
cells also function in the regime of homogeneous coverage. 

In order to discuss the motion of white blood cells in more detail, it is instructive to 
consider the parameters for a typical fiow chamber experiment. In aqueous solution and at 
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room temperature, p = g/cm^, t] = 10 ^ Pa s, and T = 2937^. Then, the dimensionless 
parameters determining cell motion become 

BAo Pe 
Pe = 4.67i?^7, / = 2.17 -, Pe, = lO.WR^Ap, At = — = 4.67i?^s, (33) 

7 7 

where R is given in pm, Ap in units of g/cm^ and the shear rate 7 in units of 1 per sec- 
onds; At is the diffusive time scale. For leukocytes in flow chambers we typically have 
i? = 5, 7 = 100 and Ap = 0.05, thus, for the two Peclet numbers we get Pe = 6 x 10^ and 
Pcz = 317, respectively. Then, / = Pcz/Pe = 0.005, that is the effect of hydrodynamic 
deterministic motion will be very strong. The experimental time scale is given by the time 
for transversing the field of view, which is about 3 s at a shear rate of 100 Hz and length of 
670 pm. The diffusive time scale At for leukocytes is about 600 s (10 min), which reflects 
their large size and shows that diffusive motion is by far not sufficient to initiate binding. 
Binding becomes more favorable in the presence of convection. For a start height of one 
radius above the wall {zq = 2), our calculations give a MFPT of about 5 s, that is much 
less than the diffusive time. However, this is still much larger than the experimental time 
scale. This proves that only those cells have a chance to bind that flow very close to the 
wall, exactly as observed experimentally. In vivo, white blood cells therefore depend also 
on other mechanisms driving them onto the substrate, including contact and hydrodynamic 
interactions with other cells. These effects have been studied in detail before. For example, 
Munn and coworkers have shown that adhesion of leukocytes close the the vessel wall in 
post-capillary venules is enhanced by red blood cells passing then>2^. King and Hammer 
have shown, using an algorithm capable of simulating several cells, that already adherent 
leukocytes can recruit other leukocytes via hydrodynamic interactions^^. The results pre- 
sented here, when specified to leukocytes, show that indeed these mechanisms are crucial 
for effective leukocyte capture under flow. 

Our results also suggest that leukocytes are sufficiently large that thermal fluctuations 
are not dominant. This changes when studying smaller particles, e.g., receptor-covered 
spheres with R ^ 1 pm, whose binding also has been investigated with flow chambers^i^i. 
Eq. ( !33|) shows that the Peclet numbers scale strongly with particle radius R, therefore, 
these beads are subject to much stronger thermal fluctuations than leukocytes. In Ref.— it 
has been verifled that indeed in equilibrium such particles obey the barometric distribution 
from Eq. f|T8l) . In Ref.— it was found that the adhesion probabihty Pad is proportional to 
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the ligand-density, pad ~ Pi- With pad ~ 1/T it follows that T ~ 1/p; as found by our 
simulations in the limit of low ligand densities. 

Throughout this paper we have considered the generic case of a constant downward acting 
force due to a density difference between the sphere and the surrounding fluid. In future work 
it might be interesting to examine also other forces which can easily be done in the framework 
presented here. As the addition formula Eq. (!26|) for falling and rotational MFPT was not 
derived under the assumption of a specific force, it is also true for non-constant forces. For 
general potential forces the falling time Eq. (!22|) has then to be replaced by Eq. (123|) . Also 
the rotational MFPT is influenced by a vertical force via the stationary height distribution. 
Neglecting gravitational force and considering only short-ranged forces like van der Waals 
or electrostatic forces would result in infinite MFPTs for the setup of the halfspace. This 
problem, however, can be solved by using an additional wall acting as an upper boundary^^. 

In this paper we assumed a rigid Brownian particle. For cells, elastic deformations might 
be relevant. For free flow, a simple scaling estimate shows that the critical value for the 
shear rate leading to substantial deviations from the spherical shape is (Eh) / (rjR)^^, where 
E = 100 Pa and h = 100 nm are Young modulus and thickness of the cellular envelope, 
respectively. The fact that the Young modulus E appears here indicates that cells tend to 
passively deform less than vesicles, whose elasticity is characterized rather by the bending 
rigidity2ii^. The scaling estimate leads to a critical shear rate of 10^ Hz, which is above 
the value of a few 10^ Hz (corresponding to Pe ~ 10^ for white blood cells) which often 
provides an upper limit in flow chamber experiments. Similar but more complicated scaling 
arguments can be made for lubrication forces which arise when the cell approaches the wall^^. 

To fully understand the rate of association between a receptor-covered particle in shear 
flow and a ligand-covered wall, our analysis should be completed by the implementation 
of an adhesion scenario, which in general should also include molecular determinants like 
residence times and receptor flexibility. If one assumes that a bond between two encountering 
molecules is formed with a certain rate, then, the MFPT for encounter as reported here 
should be a good approximation for the mean adhesion time in the limit of zero shear 
rate, because in this limit the duration of each encounter should be sufficiently long for the 
formation of an adhesion contact. Then, the proper knowledge of the MFPT could also be 
used to design a cell sorting experiment. Suppose one has a mixture of different cells each 
bearing some receptors and the wall is covered with one kind of ligand. Then, the cells are 
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flowed into the chamber and flow is stopped. Certainly, only cells that bear receptors which 
fit to the ligands can attach to the wall. If the flow is then turned on again, the attached 
cells will be separated from the other cells. If the no-flow period is much shorter than the 
MFPT, only a few cells can attach. If the no-flow period is much longer than the MFPT, 
attached cells might already start to spread and are therefore difficult to remove. Only if the 
no-flow period is of the order of MFPT one gets an appreciable number of weakly attached 
cells. In this sense our theoretical analysis might be essential for appropriate biotechnological 
applications. 
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APPENDIX A: IMPLEMENTATION OF FRICTION AND MOBILITY MATRI- 
CES 

For the numerical implementation of the friction and mobility tensors for a sphere in linear 
shear flow above a wall we use the results from Refs.— This implementation procedure 
has been described and tested in detail in Ref.— . In this appendix, we briefly summarize it 
for the sake of completeness. 

Writing the friction tensors in terms of irreducible tensors formed from 5ij, eijk, k defines 
the scalar friction functions. In the case that the normal vector to the wall is k = e^, these 
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tensors read 
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This defines the scalar friction functions (p^*^ , ip^'^ , ip^"^ , ip'''^ , (p^^ , ip^^ , (f)^'^ , ip^'^ , ip'^^ . The scalar 
friction functions (p depend only on the inverse distance of the sphere from the wall, that is 
the dimensionless variable t = R/z, which takes values from the interval [0, 1]. The friction 
functions can be expanded in powers of t. The numerically obtained first 20 coefficients of 
such a series expansion of the dimensionless scalar friction functions 

= (P'*/Q7r7]R, = ^''/67r7]R, (^^^ = (P''^/8n7]R\ 

are tabulated in Ref.^. For the other three dimensionless scalar friction functions 

the first 32 coefficients of a series expansion in powers of t are tabulated in Ref.— . For small 
values of t the series expansion converges quite well and only a few coefficients are needed 
to obtain accurate results. However, for t ^ 1, i.e., close to the wall, the friction functions 
are better described in a lubrication expansion, which reads 

4> ~ C'l^ + C2ln(l - t) + C3 + ln(l -t) + C(l - t). 

The coefficients Ci, C2, C3, C4 for the eight friction functions defined above can be found in 
Ref.—. In order to match the two limit cases, the the asymptotic expansion of the t — 1 
limit is subtracted from the friction functions 



0W = X^/nC 



ra=0 
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leading to a new series expansion: 

-Cij^^-C2 ln(l -t)- ln(l - t) 

n=l ^ \ J/ n=0 

This series is truncated at n.max = N and the coefficients Qn are calculated from the coeffi- 
cients fn, Ci. Next the coefficients Qn {n = 0, . . . ,N) are not used to calculate the Taylor 
sum, but rather to calculate the Fade approximant to this function. The Fade approximant 
is given as 



VNit) 



ao + ait + 02^^ + . . . + ttNt^ 
1 + Wt + b2t^ + ... + bNt^ 



where the coefficients the solution to 

N k 

bnON-n+k = —Qn+k-, bnQk-n = Ofc, = 1, . . . , A^. 

n=l n=l 

Finally the numerically implemented friction functions become 

m = Ci + ln(l -t)+ ln(l -t)+ VNit). (Al) 

For the calculation of the coefficients ai,bj of the Fade approximant we use the algorithm 
provided by the Numerical Recipes^^. 

Having implemented the scalar friction functions, the implementation of the mobility 
tensors proceeds by substituting ( ^ fi,(f) ^ a,ip ^ P i'i^ the above decomposition of the 
friction tensors. This defines the scalar mobility functions a**, a^^, a*^*, Z?*^*, 
Using Eq. ([H]) the dimensionless scalar mobility functions can be calculated from the scalar 
friction functions: 

= 1/0", /3" ^ 



3^ 
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In Fig. [12] we use our implementation to plot the eight dimensionless mobility functions. 
The limit of an unbounded flow corresponds to t ^ and results in 

C' = 67ir]RI, C'' = 87cr]R^I, = = = ("^ = (A2) 

where / is the unity matrix. Thus eq. ([5]) reduces to 

= 6n7]R (U - U~) , = SnrjR^ (O - n°°) . (A3) 

which are the well-known Stokes laws for the friction force and torque exerted on a sphere 
moving in a fluid with relative velocity U — U°°. For the linear shear flow considered here, 
U°° = 'jze^ and = 76^/2. 

APPENDIX B: RELATION TO THE SMOLUCHOWSKI EQUATION 

The probability distribution \1/(X, t) of a Brownian particle subject to external 
force/torque F satisfies a continuity equation dt"^ + V • J = 0. The probability flux J 
contains a diffusive and a convective part^: 

Ji = -Dijdj^H + MijFj^a (Bl) 

where D and M are diffusion and mobility matrices, respectively, and F is external force. 
In equilibrium, the flux has to vanish and the probability distribution has to become the 
Boltzmann distribution. This leads to the Einstein relation D = ksTM, which is a special 
case of the fluctuation-dissipation theorem. Using Eq. (IBip and the Einstein relation in the 
continuity equation leads to the Smoluchowski equation^S; 

dt^ = d, {M,,{kBTd,^ - F,^)) . (B2) 

We now will derive the equivalent Langevin equation. In the case of constant mobility 
{additive noise), e.g., Mjj = 6ij, the appropriate Langevin equation is given by 

dtXt = MF + gf , (B3) 

where gf is a Gaussian white noise term and the Stratonovich interpretation is used as 
explained in the main text. However if M depends on X [multiplicative noise), an additional 
drift term occurs in the Langevin equation 

a^X^ = MF + keTY + gf . (B4) 
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The following derivation of the drift term Y proceeds in two steps^. First we perform a 
coordinate transformation which makes the noise additive. In the case of additive noise the 
Langevin equation (1B3P and the Fokker-Planck equation (lB2p are equivalent. Then starting 
from the Fokker-Planck equation in the new coordinates we perform the transformation back 
to the old coordinates. Requiring the transformed Fokker-Planck equation to be of the same 
form as in Eq. (]B2|) . determines the drift term Y. 

As we use the Stratonovich interpretation for the noise process the usual rules for differ- 
entiation and integration apply and we can perform the following coordinate transformation 

X(t) 

X' = y S(X")dX", (B5) 

with some regular matrix S. The Langevin equation for the transformed coordinates then 
reads 

dtX.[ = Sdt'Xt = SMF + keTSY + Sgf . (B6) 

From the requirement that M'ij = 6ij, that is 

(SgtSgt) = 2kBTE, Eij := Sij, (B7) 

we can fix S to be the inverse of a matrix B with 

S = B-\ M = BB^ ^ Mij = BikBjk- (B8) 

As M is a symmetric positive definite matrix, it is always possible to find a matrix B with 
M = BB^. Defining 

F' := B^F + ksTSY, gf := Sgf = B^^gf , (B9) 

the new Langevin equation for the primed coordinates and with additive noise reads 

dtX.'^ = M'F' + gf . (BIO) 

The corresponding probability distribution \1/'(X', t) is the solution of the Smoluchowski 
equation 

dt<f'{-X',t) = d'MkBTd:^' - i^'^')- (Bll) 
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Next we transform (IBlip back to the unprimed coordinates. The preservation of probabihty 
requires that 



^'(X',t) = J^(X,t) 
where J is the Jacobian of the coordinate transformation^: 



J := det 



det fB) 



B 



Inserting flB12p into fIBlip gives 



(B12) 



(B13) 



(B14) 



Dividing by J we obtain for the first term on the right hand side of (]B14I) 



J-^d'kdi.m = J-\d'kd'i,J)^! + 2r^{d[J)d[^! + d'f,&^^! 



Here we made use of the identities 



j-iy'J = VB' , J-'d[J = djBj 



V = B^V 



j-'d',d'J = J-'d[{JJ-^)&J = J-\d[J)J-'&J + d[{J-^d'^J) 



kj- 



(B15) 



Again using the identity flBlSp the second term of the right hand side of flB14p can be 
evaluated to be 



J-^a^F^Jv^ = J-\diJ)F;,^ + d'.Fl,^ = a,(B,,F^v^). 
Adding both terms and inserting the definitions ( 1B8I) and ( ]B9I) we have 

dt^ = dj (kBTMjidi^ + kBTBjkidiBik)^ - MjiFi - kBTY^^) 
Comparing this with the required resuh flB2l) we can read off Y 



Y = BVB^, Y, = B,k{diB 



Ik 



Finally shifting dt^t dt^t — U°° we obtain the Langevin equation as given by Eq. (fT2|) 
combined with Eq. (fT3|) . 



28 



APPENDIX C: EULER ALGORITHM FOR A SPHERE ABOVE A WALL 



In order to solve Eq. (fT^ numerically we use an Euler algorithm. As the physical situation 
requires to use the Stratonovich interpretation of the noise term gf , the displacement AX 
of a particle from time t to time t + At depends on the position of the particle at time 
t + {l/2)At, which is not known at time t. As usual, this problem is solved by rewriting 
the Langevin equation in the Ito-version. Then the noise term can be evaluated at time t 
and as a compensation an additional drift term di{Bik)Qik is added to Eq. (fl^ ^^. Because 
B^;(9/(Bjfc) + Bikdi{Bli) = (9/(BjfcB^;) = diMu, we arrive at Eq. f|T5|) . In this equation, the 
random displacements g(At) must satisfy 

(g(At)) = 0, (g(At)g(At)) = 2MAt. (CI) 

Following Ref.^^, gi(At) is calculated from a weighted sum of normal deviate random num- 
bers Xi {xi} satisfying {xi) = 0, {xiXj) = 26ijAt. This sum is given by 



g,(At) = J]B 



where the weighting factors are the elements of the matrix B defined in flB8|) . They can 
recursively be calculated according to 

B,, = (^M^i - J2 ^ifc j ' = (^^iJ - Yl BifeB.fc j /Bjj, I > J, B^j = 0, Z < J. 

In the case of a sphere above a wall we obtain the following dimensionless weighting factors 
(cf.ii) 

Bn = \lf\ B22 = ^, 633 = v^, B42 = -B51 = (C2) 



1 



B44 = B55 = - ' = Wt|-, Bee = 1^36^. (C3) 



As pointed out in Ref.^, using the Euler method, instead of normal deviate random variables 
any uncorrelated random variable Xi — > {xi, i = 1, . . . ,6} can be chosen, as long as they fulfill 
the required relation for the first moments (xj) = 0, (xiXj) = 26ijAt. Thus, it is much faster 
to generate the random numbers according to Xi = \/12At{^i — 0.5), with = 1, ... ,6 
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being uncorrelated random variables uniformly distributed in [0, 1]. For the calculation of 
the random numbers we use the pseudo random number generator ran3 from the Numerical 
Recipes^^. 

Calculating the new configuration after each time-step using flTB]) is straightforward for 
the spatial degrees of freedom. For the update of the orientation of the sphere we use a 
coordinate system spanned by three orthonormal basis- vectors {ni|i = 1,2,3; {ni)j = 6ij}. 
The origin of this coordinate system shall be identical with the center of mass of the sphere 
and the relative orientation of this system and of the sphere are kept fixed. Given then 
an orientation update form ( |T6i) 6 := (AX4, AX5, AXg), we decompose each of the basis 
vectors into a component parallel to 6 denoted by and a component perpendicular to 
6 denoted by n_L (the index i is dropped for the sake of simplicity). These components are 
given by 

n\\ = 0{e-n), e:=e/\\e\\ 

n± = n — 9(9 ■ n). 

Then the orientation update affects only n±_ and the updated n' is given by (with 9 := \\9\\) 
n[ = 9(9 ■ ni){l - cos9) + niCos9 + 9 X nism9, i = l,2,3. (C4) 

APPENDIX D: REDUCING THE SYSTEMATIC ERROR IN MEAN FIRST 
PASSAGE TIME ALGORITHM 

Applying the Euler algorithm Eq. ( fT6l) to a mean first passage time problem gives rise 
to two sorts of errors. First there exists the statistical error, which is proportional to 
1 / V^, where N is the number of iterations the algorithm is applied. The extent of the 
statistical error of the measured mean value can be calculated during the simulation. For 
the measurements performed in sections IVIII and IVIIII typically = 10^ — 10^ iterations 
where chosen resulting in statistical errors in the range of < 1%. Error-bars in these sections 
refer to the statistical error. 

The systematic error for the mean first passage time calculated by use of an Euler algo- 
rithm scales with v^At, although the error of the particle position is only of the order of At 
— . Thus to decrease the systematic error by a factor of 10 one must increase the numerical 
cost by a factor of 100. One way to obtain accurate results at moderate numerical cost is to 
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measure the mean first passage time for various intermediate numerical time steps. Fitting 
tfiese results to a + by/At allows the extrapolation to At 0. Fig. [13] shows an exam- 
ple where this procedure was applied to the case of homogeneous coverage as considered in 
Sec. |Vl The resulting mean first passage time then deviates by 0.2% from the value obtained 
from quadrature of Eq. (l22l) . This is the same accuracy as we have for the implemented 
mobility functions themselves (cf. appendix Rl). 
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FIG. 1: Cartoon of a spherical particle with radius R moving in linear shear flow above a wall. The 
height z of the sphere center above the substrate obeys z > R. Bond formation between particle 
and wall is identified with spatial proximity between the receptor patches on the particle and the 
ligand patches on the wall being smaller than some prescribed encounter radius, that is overlap of 
the gray areas. 

FIG. 2: Falling sphere in shear flow. For different values of the shear rate (represented by the 
Peclet number Pe) and the driving force (represented by / or Pcz = f Pe) the ^-coordinate and 
the orientation angle 6 are plotted versus the x-coordinate. 

FIG. 3: Probability distribution function ^(z,t) numerically obtained from N = 10^ sample paths 
for ten consecutive points in time. The initial distribution was ^'(z,to) = ^{^ — 3) at t = to, 
Pe, = 2. 

FIG. 4: Results of first passage time simulations with encounter radius = 10^^. (a) Mean 
first passage time T as a function of Pe, for different starting heights. Dots are the results from 
simulations with = 10*^ runs and time step At = 10^^. Lines are the results from the quadrature 
of (|22p . (b) Distribution of first passage times for different values of Pe, (numerical parameters 
N = 10^ At = 10"^). 

FIG. 5: Mean first passage time dependence on the initial height zq in two dimensions. The 
sphere is covered with A^^^. = 10 receptor patches and the ligand density is pi = 0.01. We plot 
{T{zo, 9\C))0^x (+) and {T{zo, 9\C))0^x + T{z = 10|zo) (x) as a function of zq, where T{z = 10|zo) is 
obtained from Eq. (j22p . For zq > 1 + ro the latter curve is constant at the value {T{z = 10, 9\C))0^x 
as predicted by the addition theorem Eq. (j26p . (Numerical parameters: = 10^, At = 10^^.) 



33 



FIG. 6: (a) Example of a sphere restricted to move in two dimensions and covered with N,. = 4 
receptor patches, which are regularly distributed over the circumference, (b) Illustration of the 
range of 9 in which encounter occurs. This range is given by 29o with 9o{z) = arccos(2;/(l+ro))+rp. 
(c) The absorbing boundary C in the {z, ^)-plane is periodic with respect to 9 with period 9s = 
27r/iYr. For large numbers of receptor patches 9s the different patches start to overlap. Then 
encounter is possible for all values of 9. 



FIG. 7: The mean first passage time averaged over the initial orientation (log-log plots), (a) 
Plotted as a function of Pe; different symbols refer to different numbers of receptor patches, (b) 
The mean first passage time is plotted as a function of the receptor density pr oc A^^ for different 
values of Pe. (c) (r)^ as a function of A^^ in the diffusive regime (Pe ~ 0) for different values of 
the capture range ro, but fixed value of cluster-size rp = 0.001. (d) The distribution of 0-averaged 
first passage time is shown for Nr = 5, 20, 50 receptor patches. (Numerical parameters for each 
data point: N = 10^, At = 10'^. ) 



FIG. 8: (a) Illustration of the situation with a density of receptor patches pr as well as a density 
of ligands pi. The first passage time is now determined by an overlap of a receptor patch with a 
ligand patch, (b) {T)g^x as function of the Peclet number Pe and the ligand density pi for different 
values of A^^ (numerical parameters: At = 5 ■ 10~^, N = 10^). 



FIG. 9: (a) {T)g^j. is shown in the diffusion limit at Pe as a function the ligand density pi. 
Inset (plot for ~ 1): The mean first passage time scales as (T)g .j, oc 1/pf (numerical parameters: 
At = 10~^,N = 10^). (b) Dependence of {T)^^^ in the diffusive hmit at Pe on pr, pi, where pr 
has been varied by changing Nr at fixed Vp (numerical parameters: At = 10~^,N = 10^). 
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FIG. 10: (a) The mean time for a receptor to first reach a wall homogeneously covered with ligands 
{T)g was calculated as a function of the Peclet number Pe. (b) The dependence of the MFPT on 
the number of receptor patches A^^^ for different values of the capture radius tq. Lines show the 
scaling with l/Nr- (c) Dependence of {T)^^^^ on the 2D ligand density pi in the diffusive limit 
Pe ~ 0. For pi <^ 1 the mean first passage time is proportional to 1/pi (dotted lines). In the inset 
are plotted the mutual ratios of the averaged mean first passage times for A'^. = 20, 30, 70, showing 
that the dependence on the ligand density is nearly independent on the number of receptor patches 
Nj. (numerical parameters: = 10^, At = 5 • 10~^, rp = 10~^, ro = 10~^ for (a); vq = = 10~^ 
for (c)). 



FIG. 11: (a, b) Dependence of {T)^ on the receptor patch radius rp (Pe ~ 0). The dotted lines are 
fits of a/{b + rp) to the simulation results, (a) rg = 0.001, (b) rg = 0.01 (numerical parameters: 
= 1 - 3 • 10^ At = 5 • 10^5). (c) For A^^ = 30 the dependence on rp is shown for different values 
of the capture radius tq. For better comparison the ro-dependent part of the MFPT as given by 
Eq. (i22|) was subtracted. 



FIG. 12: Dimensionless scalar mobility functions. On the left the functions are plotted vs. the 
dimensionless parameter t. On the right the functions are plotted vs. 1 — t, thus better illustrating 
the asymptotic behavior for t — > 1. 



FIG. 13: The mean first passage times for Pcz = 100, zi = 1.001, zq = 2 as a function of the 
numerical time step. The points are the results from simulation experiments (error-bars denote 
their statistical error) with A^ = 10^ iterations. The full line is a fit to a + b^/At using the 
gnuplot implementation of the nonlinear least-squares (NLLS) Marquardt-Levenberg algorithm. 
Extrapolating the fit to At ^ reduced the systematic error due to the finite time step. 
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